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Abstract 

A turbulence model based on the extension 
of the algebraic eddy-viscosity formulation of 
Cebeci and Smith developed for two- 
dimensional flows over smooth and rough 
surfaces is described for iced airfoils and 
validated for computed ice shapes obtained for 

a range of total temperatures varying from 28°F 

to -15°F. The validation is made with an 
interactive boundary-layer method which 
employs a panel method to compute the 
inviscid flow and an inverse finite-difference 
boundary-layer method to compute the viscous 
flow. The interaction between inviscid and 
viscous flows is established by the use of the 
Hilbert integral. 

The calculated drag coefficients compare 
well with recent experimental data taken at the 
NASA Lewis Icing Research Tunnel (IRT) and 
show that, in general, the drag increase due to 
ice accretion can be predicted well and 
efficiently. 

1.0 Introduction 

The importance of surface roughness has 
been recognized for many years, with early 
experiments performed in relation to boundary- 
layer flows by Schlichting 1 and Nikuradse 2 . 
In these early papers, it was already evident 
that roughness might be considered with or 
without knowledge of the detailed geometry of 
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the roughness elements and that, in both cases, 
experiments were necessary to indicate the 
influences on pressure drop, skin-friction and 
heat transfer. When the elements are well 
defined, regular and with a characteristic 
dimension that approaches the thickness of the 
boundary layer, as with the ribs used to 
promote turbulence and augment heat transfer 
in some nuclear reactors, it is possible to solve 
the Navier-Stokes equations from the .wall to a 
region outside the shear layer and thereby 
provide detailed spatially-local values of 
velocity and temperature and their gradients. 
This approach is likely to require the use of 
cyclic boundary conditions, and uncertainties 
due to numerical and turbulence assumptions 
are inevitable. Where they do not meet the 
above requirements, it was not possible to 
accurately model the details of the geometry of 
the roughness elements and of the flow 
characteristics in their immediate vicinity. 

The approach to the representation of 
discrete roughness elements considered by 
Coleman and used in subsequent papers, for 
example that of Hosni, Coleman and Taylor 3 
is, in some respects, a hybrid between the 
extremes implied above. Thus, discrete 
roughness was considered in terms of a local 
element Reynolds number where the dimension 
attempted to represent the geometry. Wall-shear 
stress was then defined in two parts, one due to 
the usual form of skin friction with a blockage 



parameter and the latter based on the roughness 
Reynolds number. A similar approach was 
taken to deal with the heat transfer coefficient 
and Stanton number. These definitions, and an 
appropriate numerical description of the 
roughness parameters, were used within a 
computer code to solve numerically two- 
dimensional boundary-layer equations with 
very creditable results. Two hundred and fifty 
grid nodes were required to define the 
boundary layer at each streamwise location 
with up to 60 below the crest of the elements. 
This last feature persuaded Boyle and 
Civinskas 4 that this approach was inappropriate 
for use within their calculation method, which 
was based on the solution of thin-layer, 
Navier-Stokes equations, and it is self-evident 
that it would be less appropriate for use with 
higher-order forms of the Navier-Stokes 
equations. 

The alternative approach is to make use of 
an equivalent sand-grain roughness, as by 
Nikuradse 2 , where this quantity has been 
determined from experiment, and this was 
preferred by Boyle and Civinskas who chose to 
make use of the formulation of Cebeci and 
Smith 5 as modified by Cebeci and Chang 6 . 
This has the advantage that the turbulent flow 
calculations with surface roughness can be 
performed in the same manner as for smooth 
surfaces without the need of the 60 near-wall 
grid nodes of Hosni et al. Just as this economy 
is necessary for Navier-Stokes solvers, it is 
required with interactive boundary-layer 
methods for which the economy and efficiency 
are also important. In the particular case of the 
roughness associated with ice accretion on 
airfoils, the roughness elements cannot be 
defined and the magnitude of the effects do not 
justify attempts to follow the more expensive 
approach. 

In this paper we adopt the Cebeci-Smith 
eddy- viscosity formulation for flow over rough 
surfaces and extend it to compute turbulent 
flows over iced airfoils with a finite-difference 
method based on the solution of the boundary- 
layer equations with the interactive procedure 
described in Ref. 7. Section 2 discusses the 
method currently developed and also describes 
a new method for performing the inverse 
boundary-layer calculations which is under 


investigation. Section 3 presents results for the 
NACA 0012 airfoil with ice shapes computed 
with the fortified LEWICE code 8 and the total 
drag coefficients computed by the interactive 
boundary-layer (IBL) method developed by 
Cebeci 7 and described in Refs. 8 and 9. The 
paper ends with a summary of the more 
important conclusions. 

2.0 Turbulence Model for Iced Airfoils 

We use the calculation method of Ref. 7, 
which is based on the solution of the two- 
dimensional incompressible boundary-layer 
equations, which employ the eddy-viscosity 
concept to model the Reynolds shear-stress 
term 


du dv 

s + ^ =0 


(hi chi du e d ,'du . 


( 1 ) 

( 2 ) 


where b=\ + e m /v . On the airfoil, Eqs. (1) 
and (2) are subject to the boundary conditions 

y=0, m=v= 0 (3a) 

y =8, u =u e (x) = u° (. x ) + 8u e (x) ( 3 b) 


where u°(x ) denotes the inviscid velocity and 

8u e (x) the perturbation velocity due to 
viscous effects which is given by 


it J x a da x— a 


( 4 ) 


in the interaction region ( x a , Xb). For the 
calculation of wake flow, Eqs. (1) and (2) are 
subject to conditions on the dividing centerline 

y c 


y=y c ,v= 0 (5a) 

and to freestream conditions 
y — > ± °°, u -+u e (x) = u° e (x)+8u e (x ) (5b) 
As discussed in Ref. 10, Eqs. (3b) and (4) 
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can be written as 


n 

u e (x) = u° e (x)+'£c ij ( u e d*)j 

;=! (6) 

Here c ij denotes the interaction coefficient 
matrix, which is obtained from a discrete 
approximation to the Hilbert integral in Eq. (4). 
In this form, Eq. (6) provides an outer 
boundary condition for the viscous-flow 
calculation, which represents the viscous- 
inviscid interaction. It can be generalized to the 
form 

= (*)+£<*/ [(w« 8* )j — (m, fi* 

' =1 (7) 


where w = u e /u a > and, with c u denoting the 
interaction coefficients, 



In the above equations, / denotes a 
dimensionless stream function defined by 

¥= f (<%, Tj) ( 12a ) 

and primes denote differentiation with respect 
to Tj given by 

V =JiUvxy ( 12 b) 


where u?(x) corresponds to the inviscid 
velocity distribution which contains the 

displacement thickness effect (<5*) x ‘ computed 
from a previous sweep, as discussed in Ref. 
10 . 


The numerical solution of the above 
equations for a given eddy-viscosity 
formulation was obtained by using Keller's 
box scheme after the continuity and momentum 
equations and their boundary conditions are 
expressed in the following form 


In Eq. (9b) the parameter gi results from the 
discrete approximation to the Hilbert integral 
and is given by 


«-i 

8i = ™ K + X °ij C Dj -Df)-c u D? 

(13) 

where 


D = 



( VeW-fe ) 


(14) 





( 8 ) 


On the airfoil 


*7 = 0 , /=/' = 0 


(9a) 


The expression for gi on the wake is nearly 
identical to that for the airfoil, Eq. (13), except 
that now Eq. (14) is given by 

D = <— )> [h- (J). - IJ -. ) - if, )] 

J (15) 


7 1= Ve> W ~ Cu (7 leW~f e ) = gi (9b) 

and in the wake 
V=l-e, f=w 
77=0, /= 0 

77 = 77 ,, f-w ( 10 ) 

W - Cu [w (77, -7 l-e)-(f e -/_, )] = gi 


In regions of flow separation, the so-called 
FLARE approximation is introduced in which 

the convective term, / (df l d^) is set equal to 
zero in the recirculating region. A nonlinear 
system of algebraic equations resulting from 
the finite-difference approximations is 
linearized by Newton's method and solved by a 
block-elimination procedure. The details of the 
solution procedure are available in several 
references, for example Ref. 1 1. 

For flow over a smooth surface, the 
algebraic eddy-viscosity formulation of Cebeci 
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and Smith is given by 


f (£m)i=L 2 


du 

dy 


Yir 


£m ~ 


<>($»)* = 0.0168 m , Fy, 


tr 


(16a) 

(16b) 


Ay = 


.0.9 (v/m t )[/^-^ + exp(-#/6)] 

5 <kt< 70 

( 20 ) 

' 0.7(v/u T )(#)°- 58 70 < kt < 2000 


where L denotes a modified mixing length 

L = 0.4y [1 - exp(-y/A)] (1?) 

and A is a damping-length parameter 

A = 26vuj\ u x = (T/p)i/« 

The parameter y, r represents the transition 
region 


where « T represents the friction velocity given 
by u T = J%Jp . 

The Cebeci-Smith eddy-viscosity 
formulation with L expressed as Eq. (19) has 
been applied for flows over iced airfoils by 

providing a link between k£ and the 
roughness of the airfoil surface with ice. This 
method described below was used for the 
calculations reported in this paper. 


y tr = l-exp[-G (x-x tr ) f* — 


(18a) 


where x tr is the location of the beginning of 
transition and G is defined by 




-1/34 

*/r 


(18b) 


with Rx, r — (UeX/ v)tr and C is a constant equal 
to 60. 


The Cebeci-Smith eddy-viscosity 
formulation has been extended for flow over a 
rough surface (Ref. 6) by expressing the 
modified mixing length L as 


The integral boundary-layer method used to 
determine the heat transfer coefficient in the 
LEWICE code (Ref. 12) accounts for the 
surface roughness of an iced airfoil by 
expressing the equivalent sand-grain roughness. 
k s as a function of liquid water content (LWC), 
static air temperature (T s ), and airspeed ( ). 

With c denoting the airfoil chord and for 
( k s /c)base = 0.001 17, the equivalent sand-grain 
roughness k s is expressed in the following 
form 


jr — r -i f k s /c 1 

K s — l~ ; r iLwC'i — — ; — : Jr, 


( k s / C )ba 


{ k s /c ) b a 


r kgfc 1 k s . 

•1/ • < \ JV--( )base-C 

\k s /c) base c (21) 


L = 0.4(y + Ay){l-exp[-(y + Ay)/A)]} 

(19) 

recognizing that the velocity profiles for 
smooth and rough walls can be similar, 
provided that the coordinates are displaced. In 

Eq. (19), Ay is expressed as a function of an 
equivalent sand-grain roughness 

parameter £/(=^m t /v) , i. e .. 


where each sand-grain roughness parameter is 
given by 


^ Uc 

( k s /c )base 


]lwc = 0.5714 + 0.2457 (LWC) 


+ 1.2571 (LWC) 2 


j- ks/c 
( K/C )base 


= 0.047 7; - 11.27 


(22a) 

(22b) 
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I tttt ]v- = 0.4286 + 0.0044139 V. 

( k s /C )basc 

(22c) 

These expressions, good for both rime and 
glaze ice, are empirical, are based on 
experimental data reported in Ref. 12, and do 
not account for the effect of time on the ice 
roughness. The experimental data of Olsen et 
al. 13 shows that, for glaze ice condition, the 
roughness increases with time, rapidly at first 
then more slowly. Their data also shows that 
rime ice is never as rough as glaze ice. 


Recent studies conducted by Shin et al. 9 

show that k s is independent of airspeed 
but is also a function of the median volume 
droplet (MVD) size as well as a function of the 
parameters in Eq. (21). As a result, they write 
the expression for sand-grain roughness, given 
by Eq. (13) as 


K = 0.6839 [ 


k s /c 
( k s !c )ba 


•L^c.r 


k jc 

( k s / C )ba 


-]7i 


J- ks/c 

( k s jc )base 


)base-C 

C 


(23) 


with 


r k S /C 
(ks/c) b 

ase 


]mvd - 


1 MVD < 20 

1.667-0.0333 MVD 
MVD > 20 
(24) 


new method which not only can make the 
boundary-layer calculations more efficient but 
can also lead to substantial savings when used 
in the Navier-Stokes calculations is being 
considered and the approach is described 
briefly here. This new method also uses the 
Cebeci-Smith eddy-viscosity formulation but 
replaces the true boundary conditions at y = 0 
by new "wall" boundary conditions defined at 
some distance y 0 outside the viscous sublayer 
to avoid integrating the equations through the 
region of large y gradients near the surface. 

We choose this y 0 to be the distance 
given by 


yo = (v/u r )y + 0 (25) 


with y 0 being taken as a specified constant. 
In that case, the "wall" boundary conditions for 
w at y = y Q can be represented by the law-of- 
the-wall expression given by Eq. (6.34) in 
Cebeci and Bradshaw 14 

u a = u-c [-In (y„ + Ay + ) + c-Aw + ] 

Lx* J (26) 

For Prandtl-Schlichting sand-grain roughness 
(see Cebeci and Bradshaw 14 ) the shift in the 
velocity profile, A u + , can be represented by 


Au + = 


,3.94 - 6.276 In# + 2. 876 (In#) 2 
-0.290 (In#) 3 4.4 <#<70 


l- In# -3.70 

K 


#>70 


The present method of applying the Cebeci- 
Smith eddy-viscosity formulation for an iced 

airfoil is to use Eq. (16) with L and A y given 
by Eqs. (19) and (20) and with k s given by a 
modified form of Eq. (23). With the eddy- 
viscosity formula defined in this way, Eq. (8) 
can be solved subject to its boundary 
conditions given by Eqs. (9) and (10) in an 
.. inverse mode. 


(27) 

It can then be shown from the continuity 
equation that the second "wall" boundary 
condition is 


du Tr u 0 * . 
v* =-y<>-r L {—+-( 
dx u T 


-(— ) In [1 + — - — ] 
Kyo (A y/y o y 


Up to this point, the current turbulence 
model for an iced airfoil has been discussed. A 


.+ d(Au+) 

1 (28) 
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where 


where 


-6.276 + 5.753 In# 
-0.870 (In#) 2 


, + d(Au+)_ f 
5 dkt 

1 

^ — 
K 


4.4 < # <: 70 


# > 70 (29 ) 


The changes in shear stress between y = 0 
and y =y Q can be computed from 



where 


(30) 


Jul 


du. 


a = 0.3— — u e — - — 


dx 


dx 


(31) 


In terms of transformed variables defined 
by Eq. (12), the "wall" boundary conditions at 

V = Vo can be written as 

fo=—{~ In [ (— )//?T (Vo + Ari 0 ) ] 

Woo ^ Wo© 

+c - A#} (32) 

_j_ fo^ 1 n ^ Vo du z j. j Moo ^ jf 

og 2 2 w« dx w T 


+ 1 At? 

*Vo 


In (1 + 


1 

A 77/77 * 



dAu + 

dkt 


(33) 


Similarly, Eq. (30), which provides the 
relationship between the wall shear and that at 
y-y<y> can be written as 


2 w 2 r , xu e du e 

= — {#& + 77* — —A 
ul dx 




-o.s^^-^) 2 -^/;/;]} 


Vo = v^(— ), 

* 



(35) 


We note that when the momentum 
equation, Eq. (8), is solved subject to the 
"wall" boundary conditions at y = y a , namely, 
those given by Eqs. (32) and (33), the 
modified mixing length expression given by 
Eq. (19) is now 


L = 0.4 (y + Ay) (36) 

The resulting inner eddy-viscosity formula in 
Eq. (16a), expressed in transformed variables 
without the intermittency term, is 

te* )t = — = 0. 16(77. + A7j) 2 k I JR X 

v II (37) 

The solution of Eq. (8) subject to those 
boundary conditions given by Eqs. (32) and 
(33) is somewhat more involved than the 
solution with "true" wall boundary conditions 
since the "new" wall boundary conditions are 
nonlinear and friction velocity is not known. 
An efficient and accurate solution procedure 
requires not only the linearization of the 
boundary conditions with Newton's method as 
employed in Keller's box scheme, but also the 
coupling of the solutions to the wall shear 
given by Eq. (34). 


A convenient numerical approach for the 
solution procedure is to treat the friction 
velocity as an eigenvalue and obtain the 
solutions of Eq. (8) subject to Eqs. (32) and 
(33) for successive values of m t . More 
specifically, if we let m t = w, then the 
successive values of w can be determined from 
the following expression 


w 


v+l 


w v - 


ft(w v ) 

( d<p!dw ) v ’ 


v=0, 1 , 2 ,... 

(38) 


(34) 
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based on Newton's method. In the above 

equation, 0 represents essentially the error term 
resulting from Eq. (34), that is 



0(W V ) = (w v ) 2 - [RHS of Eq.(28)] = 0 ( 39 ) 

The term dfy/dw is obtained by solving a set 
of variational equations obtained by 
differentiating the difference equations of Eqs. 
(8), (32) and (33) with respect to w. 

A better and perhaps a more efficient 
numerical approach for the solution procedure 
is to treat w as an unknown and noting that 


2*1 (40) 

solve the boundary-layer equations as a fourth- 
order system rather than a third-order system. 
Both numerical approaches are under 
investigation. 

3.0 Results and Discussion 

The current method of applying the Cebeci- 
Smith eddy-viscosity formulation described in 
the previous section was used to make 
calculations for the same icing conditions used 
in the recent test in the Icing Research Tunnel 
(IRT) of the NASA Lewis Research Center 15 , 
and comparisons were made in predicting the 
measured ice shapes on the leading edge of an 
NACA 0012 airfoil together with their total 
drag coefficients. Before we present a 
comparison between measured and calculated 
results, it is useful to present a brief summary 
of the experimental data reported in more detail 
in Ref. 15. 

The measurements were made in the NASA 
Lewis Icing Research Tunnel on a 6 ft span, 21 
in. chord NACA 0012 airfoil with a fiberglass 
skin. The model was mounted vertically in the 
center of the test section as shown in Fig. 1. 
During all icing runs, the model was set at an 

angle of attack of 4°. The icing conditions 
reported in Ref. 15 can be grouped as (1) low 

speed with V, =150 mph and a high liquid 
water content, LWC of 1.0 g/m 3 , and (2) high 
speed, V« = 230 mph and low LWC = 
0.55 g/m 3 . Water droplet size, MVD, was kept 

constant at 20 |im in both groups. The ice 
accretion time in group one was 6 minutes and 
in group two was 7 minutes. Total 


temperatures in both groups were varied from 

28°F to -15°F to cover glaze, rime and 
transition regimes. 


3.1 Comparison Between 
Measured Ice Shanes 


Calculated and 


Calculations were made with the fortified 
LEWICE code for the comparisons with the 
experimental data. The time step used in the 
calculations for ice accretion was one minute. 


Figure 2 shows a comparison between the 
calculated and measured ice shapes at an air 
speed of 150 mph for seven total temperatures 
in which the ice shape changes from white, 
opaque rime ice to slushy, clear glaze ice when 
the temperature is increased. At low 
temperatures at which the ice shapes 
correspond to rime ice, the comparison 
between calculated and measured ice shapes are 
good, consistent with the previous study 
reported in Ref. 9. Icing limits are predicted 

well for temperatures below 18°F. At warmer 
temperatures the calculated results predict more 
run back, which results in more ice accretion 
beyond the measured icing limits. While the 
calculated results predict the direction of hom 
growth of the glaze ice, in general the size of 
the predicted ice shape is larger than the 
measured size. 


Figure 3 shows a similar comparison for 
the airspeed of 230 mph. As in Fig. 2, 
computed results agree well with measurements 

at all temperatures, except at 28°F where an 
overprediction of upper hom ice is seen. 


3.2 Comparison Between Calculated and 
Measured Drag Coefficients 


A comparison between calculated and 
measured drag coefficients can be made with 
the interactive boundary-layer method of Refs. 
8 and 9 using the turbulence model given in 
this paper for an ice shape either given by 
experiment or computed by the fortified 
LEWICE code. The former choice is usually 
not practical for two main reasons. 
Experimental ice shapes contain discontinuities, 
and the flowfield calculations require smooth or 
near smooth shapes. The second, perhaps the 
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most important reason, is the requirement that 
the extent of ice on the airfoil be known and 
this is difficult to measure or specify. For 
these reasons, the calculations here were made 
using the ice shapes determined by the fortified 
LEWICE code, as described in the previous 
subsection. 

A previous study 9 , employing an approach 
similar to the study conducted for the Olsen et 
al. data 13 taken on a NACA 0012 airfoil, 
showed that while the calculated drag 
coefficients predicted the correct trend with 

change in total temperatures ranging from 32°F 

to -15°F, the magnitude of the drag coefficients 
were lower than the measured ones, especially 
for the rime ice. This difference may be due to 
(1) the turbulence model and/or (2) the 
roughness associated with the iced surface. 

While the present turbulence model has 
been evaluated and found to be satisfactory for 
turbulent flows over rough surfaces, its 
accuracy still remains to be explored for iced 
airfoils, since this model primarily uses the 
correlations developed for high Reynolds 
number flows over flat plates with roughness. 
The correlation between the drag of a flat plate 
covered with uniform sand-grain roughness 
and the drag of irregular ice shapes on airfoils 
is not at all well established. The ice attached to 
the leading edge is in a flowfield in which the 
inviscid surface velocity varies between the 
extremes of zero at the stagnation point and the 
maximum velocity, and ice particles may be 
subject to what is known as the drag 
magnification. This concept was advanced by 
Nash 16 and it was confirmed experimentally 
that the drag of objects placed in a velocity field 
is not simply proportional to the local dynamic 
pressure but varies according to a different law, 
dependent on the initial momentum loss and 
some power of the local velocity. Thus the 
drag contribution of roughness placed in a high 
velocity region may be more than that deduced 
from the application of flat-plate values locally. 
Not accounting for the drag magnification in 
the turbulence model may lead to drag values 
lower than those observed experimentally. 

A numerical study was conducted to 
investigate ways to compensate for this 
deficiency in the roughness model (i.e. drag 


magnification). The current IBL method 
calculates the drag coefficients with the 
roughness parameter (k s )mu given by 

(kr)lBL = 2(kr)E q .(23) (41) 

applying it only over the surface of the ice. For 
this numerical study, the code was modified to 
allow for roughness on both the ice and the 
airfoil surface downstream of the ice to 
investigate the effect of the extent of roughness 
on drag. The results showed that agreement 
between calculated and measured drag 
coefficients for rime ice shapes became much 
better by extending the range of the roughness 
on the airfoil surface and placing a lower limit 
of kjc = 0.002 on the equivalent sand grain 
roughness, which otherwise would become 
very small for rime ice. The extent of the 
roughness which resulted in the best agreement 
with the experimental drag coefficients for rime 
ice shapes was found to be 50 percent of the 
airfoil chord, and this extent was used in all 
drag calculations presented in this paper. 

Table 1 and Figs. 4 and 5 show a 
comparison between the calculated and 
measured drag coefficients for the ice shapes 
shown in Figs. 2 and 3, respectively. As can 
be seen, the values of the drag coefficient 

remain nearly constant up to around 12°F for 
both speeds before they increase sharply with 
increasing temperature as the ice shape changes 
to that of glaze ice. The adjustment made by 
extending the range of roughness on the airfoil 
allows the calculated drag values to agree well 
with experimental data at low temperatures. 
The computed drag coefficients begin to rise 

sharply at around 18°F and reach a peak at 

22°F for Yo = 150 mph. Up to this 
temperature, the trend of calculated results and 
their magnitudes are in good agreement with 
measurements. However, for temperatures 
higher than 22°F, while the computed drag 
coefficients begin to decrease, the measured 

ones continue to rise sharply until 28°F before 
they start decreasing rather sharply. For \C = 
230 mph, however, the calculated results does 
a good job of following the trend in measured 
values. A typical computing time (CPU time 
only) for a calculation with multiple time steps 
was less than 50 seconds on a CRAY X-MP. 
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Figure 6 shows the variation of the drag 
coefficient of the airfoil with angle of attack for 

a fixed ice shape computed at a = 4°, V, = 
150 mph, T t = 22°F (Fig. 2c). The calculated 
drag coefficients for angles of attack from 0° to 

4° agree very well with experimental data, and 
the calculations exhibit no numerical 

difficulties. In the angle of attack range of 4° 
to 8° it was only possible to compute drag 
coefficients for the values of a indicated in Fig. 

6. For intermediate values of a, the numerical 
solutions for wake flow did not converge 
properly, so the calculations could not be 
completed. The cause for this behavior of the 
code is under investigation. Except for the 
difficulties in convergence at these intermediate 

values of a, the drag coefficients computed 
without any convergence problems and shown 
in Fig. 6 indicate a very good agreement with 
measured values. 

4.0 Concluding Remarks 

The results presented in this paper, as those 
in a previous paper (Ref. 9) indicate that a 
combination of LEWICE and IBL provide an 
efficient and reasonably accurate tool for 
predicting ice shapes on airfoils and for 
determining their effect on the increase of drag. 
Improvements are needed, however, in the 
interactive boundary-layer method so that the 
occasional lack of convergence of the solutions 
at angles of attack near and/or post-stall will be 
resolved. Studies are also needed to better 
estimate the extent of roughness on the airfoil 
surface. Finally, further work is needed to 
determine why predicted ice profiles are larger 
than the measured profiles at the warmer 
temperatures. 
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Table 1 . — Effect of total air temperature on drag coefficient, 
(a = 4°, MVD = 20 ^m) 


Total 

temperature, 

°F 

Experimental 

drag 

coefficient 

Calculated 

drag 

coefficient 

28 

0.0578 

0.0346 

25 

.0540 

.0372 

22 

.0315 

.0392 

18 

.0271 

.0351 

12 

.0229 

.0217 

1 

.0229 

.0209 

-15 

.0233 

.0202 


(a) V„ = 1 50 mph, LWC = 1 .Og/m 3 , t = 6 min. 


Total 

temperature, 

°F 

Experimental 

drag 

coefficient 

Calculated 

drag 

coefficient 

28 

0.0428 

0.0470 

25 

.0371 

.0294 

22 

.0311 

.0202 

18 

.0268 

.01 95 

12 

.0255 

.01 95 

1 

.0234 

.0195 

-15 

.0218 

.0192 


(b) V„ = 230 mph, LWC = 0.55g/m 3 , t = 7 min. 
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Figure 2. — Effect of total air temperature on ice shape for the experimental conditions in Table 1 (a). 
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(a)T t =28«F 
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Figure 3.— Effect of total air temperature on ice shape for the experimental conditions in Table 1 (b). 
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Figure 4. — Effect of total air temperature on the airfoil drag 
coefficient for the experimental conditions in Table 1(a). 
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Figure 5.-Effect of total air temperature on the airfoil drag 
coefficient for the experimental conditions in table 1(b). 
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Figure 6. — Variation of the drag coefficient of the airfoil with 
angle of attack with ice shape corresponding to glaze ice 
shown in Figure 2 (c). 
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